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Abstract — Nonequilibrium quantum mechanics can be solved 
with the Keldysh formalism, which evolves the quantum mechani- 
cal states forward in time in the presence of a time-dependent field, 
and then evolves them backward in time, undoing the effect of the 
time-dependent field. The Feynman path integral over the Keldysh 
contour is employed to calculate the strongly correlated Green's 
function. We examine the accuracy of this procedure for the sim- 
plest problem that requires a nonequilibrium formulation: the f- 
electron spectral function of the spinless Falicov-Kimball model. 



I. Introduction 

ELECTRONS are correlated when the Coulomb repul- 
sion between them is strong enough that it plays a 
significant role in determining the motion of the elec- 
trons through the crystal. Correlated electrons are of 
interest to the military, because their properties (metal- 
lic/insulating/magnetic/superconducting/etc.) can be easily 
tuned by changing pressure, chemical composition, irradiation 
with electromagnetic fields, and so on, and form the basis of 
many so-called smart materials and devices. In addition, a large 
number of materials of interest to the military (like Plutonium) 
have strong electron correlations 

Our research problem involves understanding strongly cor- 
related materials when they are placed under intense electro- 
magnetic fields that can drive them out of equilibrium creat- 
ing interesting dynamical and relaxational effects (examples in- 
clude intense pulsed laser irradiation or interacting with large 
amplitude microwaves). Our main focus is to the Josephson 
junction device (0 (a sandwich of a superconductor-barrier- 
superconductor which has the potential for ultra high speed 
digital electronics |2|), but the principles can be applied to a 
large number of different devices. The many-body formalism 
for nonequilibrium problems is solved by a Feynman path inte- 
gral over the so-called Keldysh contour 1 3 1, 1 4 1 which involves 
an evolution forward in time as the external fields are turned 
on, evolution out to a long time, then an inverse evolution back- 
ward in time where the fields are turned off. Solving the Feyn- 
man path integral requires evaluating finite-sized determinants 
of discretized matrices that represent the continuous matrix op- 
erator along the contour. We typically require the determinant 
of approximately 500 general complex matrices with sizes up 
to about 2100 x 2100 for a production run. This computa- 
tional effort is easily parallelized. The numerical solutions suf- 
fer from a discretization error that gets worse as the temperature 
is lowered, and accurate calculations require a careful extrapo- 
lation with different discretization sizes to the limit where the 



discretization goes to zero. Here we benchmark the numerics 
by solving for the /-electron spectral function of the spinless 
Falicov-Kimball model |5J- 

The Falicov-Kimball model is the simplest model of elec- 
tron correlations (and the problem we investigate is the sim- 
plest nontrivial Keldysh problem). It possesses two types of 
electrons: conduction electrons, which can hop to any of their 
nearest neighbors and localized (/) electrons which are local- 
ized on the lattice sites. There is a Coulomb repulsion U be- 
tween conduction electrons and localized electrons on the same 
lattice site. As U increases, the conduction band splits into two 
bands with an energy gap in between, and undergoes the so- 
called Mott metal-insulator transition. In this contribution, we 
restrict ourselves and all formulas to the case of half filling, 
where half of the lattice sites are occupied by the conduction 
electrons and half by the localized electrons. In this case, the 
Hamiltonian becomes [5 1 

(»i> 1 

- |E( c ^+/M), (i) 

i 

where c\ (Cj) creates (destroys) a conduction electron at site 
i, fj (fi) creates (destroys) a localized electron at site i, U is 
the interaction strength, and t* is the hopping integral. The 
symbol Z represents the number of nearest neighbors, and (ij) 
denotes a sum over all nearest neighbor pairs. The first term is 
the kinetic energy of the conduction electrons, the second term 
is the Coulomb repulsion, and the third term is the chemical 
potential times the filling for both electrons. 

We solve this problem on an infinite coordination number 1 6 1 
(Z — > oo) Bethe lattice, which has a noninteracting density of 
states (DOS) that is a semicircle 

P(e) = ^V4** 2 -e2. (2) 

The noninteracting bandwidth is At* . We choose t* as our en- 
ergy unit and set it equal to 1 . See Ref. [71 for a review of the 
equilibrium and linear response solutions. 

II. Formalism 

We start with an examination of the retarded local Green's 
function for the conduction electrons, defined to be 

G n (t) = -z9(t)Tv(e-^{c 3 (t),c](0)} + )/Z, (3) 



with 8(t) the unit step function, Tr denoting the trace over all 
conduction electron and localized electron states of the lattice, 
j3 = 1/T, Cj(t) = exp(ifH)cj exp(— itH), the braces denote 
an anticommutator, and Z = Tr(exp(— /3H)). We chose to 
evaluate the local Green's function at site j, but it is the same at 
every site when there is no long range order. In the limit of in- 
finite coordination number, we find that the many-body self en- 
ergy for the Green's function becomes local |6|, and the many- 
body problem for the lattice can be mapped onto the many -body 
problem for an impurity in a time-dependent field with a self- 
consistency relation to the lattice |8|. The Fourier transform of 
the retarded Green's function can be determined (on the Bethe 
lattice) by solving a simple cubic equation |9], 1101 
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)G(w) - u = 0, (4) 



where one must choose the (causal) physical solution deter- 
mined by the root with a negative imaginary part (when the 
imaginary part is nonzero) and by continuity when real. The 
conduction electron DOS is defined by A(uj) = — lmG(aj)/iv, 
the electronic self energy (on the Bethe lattice) satisfies 
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and the dynamical mean field X(ui) is defined to satisfy 



G(w) = 



1 
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(5) 
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which can be thought of as the Fourier transform of the time- 
dependent field for the impurity problem (indeed the dynamical 
mean-field theory approach is to construct an impurity problem 
in a time-dependent A field and then adjust the field until the 
Green's functions for the impurity are equal to the local Green's 
functions for the lattice J8|). 

In order to calculate the /-electron Green's function, we must 
first start with the impurity problem, whose Hamiltonian is the 
same as the lattice Hamiltonian, but there is no hopping term, 
since it is restricted to the single site of the impurity. The hop- 
ping of the conduction electrons is mimicked by the time de- 
pendent A field which destroys a conduction electron at time 
and creates a conduction electron at time t with "strength" A(t). 
The (greater) /-electron Green's function is then defined by 

G>(t) - -Tr(e-^ n ^S c (\)f(t)p(0))/Z vmp , (7) 

with f(t) = exp(itHim P )f exp(—itTi.im P ) and the evolution 
operator given by 



S C (X) = T c exp 



dt / dt'J(t)X c {t,P)c{P) 



(8) 



The time-ordering is along the Keldysh contour (see Fig. 
and the contour-ordered dynamical mean field is found from a 
Fourier transform of X(oj) 
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Fig. 1 

Keldysh contour for evaluating the /-electron Green's 

FUNCTION AT TIME t. THE CONTOUR RUNS FROM t = TO t = t, THEN 
BACK FROM t = t TO t = AND FINALLY GOES ALONG THE IMAGINARY 
AXIS DOWN TO t = — i/3. THE HELD A c - Xt IS ACTIVE ON THE FORWARD 
BRANCH, AND THE FIELD A c IS ACTIVE OVER THE BACKWARD BRANCH 
AND THE IMAGINARY BRANCH. WHEN WE DISCRETIZE THE MATRIX 
OPERATOR OVER THE KELDYSH CONTOUR, WE EVALUATE THE 
INTEGRALS VIA A RECTANGULAR (MIDPOINT) SUMMATION WITH A STEP 
SIZE OF At REAL ON THE REAL AXIS AND At IMAG ON THE IMAGINARY 
AXIS. WE CHOOSE AtiMAG = 0.05 AND VARY At REAL FROM 0. 1 TO 
0.0125 IN OUR CALCULATIONS. WE TYPICALLY USE NO MORE THAN 
1000 TIME STEPS ON THE OUTWARD BRANCH OF THE CONTOUR. 



where /fzj(w) = 1/[1 + exp(/3u>)] is the Fermi-Dirac distribu- 
tion and 9 c (t — i') = if i' is in front of t on the contour c and 
1 if it is behind. 

This Green's function can be solved directly by evaluating a 
Feynman path integral over the Keldysh contour to yield II II . 
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S c {i-t')+ / dt"g aux (t,t")X(i",F) 



which involves the determinant of a continuous matrix opera- 
tor (note the path integral is the time-ordered product along the 
integration contour, which is the Keldysh contour here). The 
function g aux appears in Table U It is the Green's function for 
a noninteracting Fermion that evolves in a time-dependent field 
Xt(i, P) — —iU8 c (t — P)6 c (t — P), which arises from the time 
dependence of the localized electron II II . J7). In this sense, 
one has the fields A c — Xt acting on the forward branch of the 
Keldysh contour, and the field A c acting on the backward and 
imaginary branches of the contour. The (equilibrium) greater 
Green's function satisfies a spectral formula with the DOS and 
the Fermi-Dirac distribution function (because the equilibrium 
distribution function is known) 



G>(t) = 
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[f FD (u)-l]A f (u). 



(11) 



Using the fact that the greater Green's function satisfies 
G>(i) = G^*(-t) and the fact that the DOS at half filling 
is an even function of li (due to particle-hole symmetry) yields 
the final equation for the /-electron DOS 111 II . Q 



dtRe{G^(t)}cos(ujt). (12) 



TABLE I 

g aux (t, P) FOR DIFFERENT ORDERINGS OF t, t, AND ¥ ALONG THE 

contour c. The symbol f satisfies fo = 1/[1 + exp(iUt - /3U/2)]. 



Green's function value 



domain 



6,exp[itf(t-f)/2] 
£pexp{iUt-iU(t + t')/2] 
iaexp[-iU(t-F)/2] 
(&-l)exp[itf(i-F)/2] 



t <t <t' 
t<t<i' 
i<¥ <t 
t <i' <i 



(Co ~ 1) exp[-i£7i + if/(t + f')/2] t' <t<t 
(go - l)exp[-zL7(f- F)/2] i' <t<t 

In addition to computing the Green's function on the real fre- 
quency axis, one can also compute it on the imaginary (Mat- 
subara) frequency axis, at the Matsubara frequencies ioj n = 
iirT(2n +1). In this case, one can formulate the expressions as 
the determinant of a discrete matrix, so there is no error asso- 
ciated with discretizing a continuous matrix operator 11 II . 171 . 
Since the Matsubara Green's function can also be expressed as 
an integral over the DOS 
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(13) 



we have an independent way to verify the accuracy of the DOS 
by comparing the integral formula for G(iu) n ) with the result 
directly calculated on the imaginary axis. Usually the accuracy 
is worst for the lowest Matsubara frequency. 

There also are a number of moments and properties of the 
DOS that can be tested. First, the DOS is always nonnegative — 
negative values of the DOS for some region of frequency indi- 
cate an error in the calculations. Second, one can work out 
explicit values for the first three moments. At half filling, these 
satisfy 



d,LoAf(uj) = 1, 



(14) 



duA f (u)uf FD (u) = LT(( c t c /t/) - _), (15) 
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In practice we add a small shift (always less than 0.006 in mag- 
nitude) to the DOS in order to satisfy the unit weight sum rule 
[Eq. 1141 1. Then we check the accuracy of the other two sum 
rules [by independently calculating the correlation function on 
the right hand side of Eq. dl5H . At t = 0, we can use the mo- 
ments in Eq. ( 1 1 414 1 6I > to show that 
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(17) 
(18) 
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Thus the first derivative depends on temperature, but is purely 
imaginary, while the second derivative is real and independent 
of T. Hence we expect ReG^* (t) to depend weakly on T for 
small times, and to show stronger dependence at large times. 



III. Computational Algorithm 

The main effort of this nonequilibrium calculation is to com- 
pute the /-electron Green's function, which requires the deter- 
minant of the continuous matrix operator in Eq. dlOi . To cal- 
culate this, we must first decide on a discretization along the 
contour to evaluate line integrals over the contour (see Fig.Q. 
Our choice is to use a step size of Ai rea j on the real axis and 
A<j ma g = 0.05 along the imaginary axis. We use a midpoint 
rectangular integration rule, evaluating the function at the mid- 
point of each discretized piece of the contour, for the approxi- 
mation to the line integral. Hence 



f dih(t) = Y l W j h(i j ), 

J c 



(20) 



with the weight function Wj satisfying Wj — ±Ai rea j on 
the real axis (positive on the forward branch and negative on 
the backward branch) and Wj — — *Aij ma g on the imag- 
inary axis. The times where the function is evaluated are 
tj = (j — l/2)At rea j on the forward branch, tj — (2jmax 
/ + l/2)At rea j on the backward branch, and tj = i(2jmax — 
j + l/2)Aij ma g on the imaginary branch (jmax is the number 
of points on the forward branch of the Keldysh contour). Using 
this scheme, a Dirac delta function is approximated by 



S c (tj 



(21) 



The evaluation of a discrete approximation to the determi- 
nant of the continuous operator is now a straightforward proce- 
dure 1 1 1 1 . First we note that 

Det c (l + M) = cxp(Tr c [ln{l + M}]) = cxp(^ -Tr c M"), 

n 

(22) 

is a relation relating the determinant to the trace of a series of 
powers of the matrix M. The symbol Tr c denotes the trace of a 
matrix operator over the Keldysh contour, and it satisfies 

Tr c M = / dtM(i, t)=Y^ WiM(ti,U). (23) 

Hence the trace of the powers of M becomes 

Tr c Af"= W il ...W in M(t il ,t i2 )...M(t in ,t il ). (24) 



Now we define a new discrete matrix to satisfy M(ti,tj) = 
WiM(ti, tj). Then we find the trace in Eq. ( 1241 becomes 

Tr c M"= M(t il ,t i2 )M(t i2 ,t i3 )...M(t in ,t il ) = TrM n , 

ii ■ ■■i n 

(25) 



and leads to the final formula for the determinant 
Det c (l + M) =Det(l+M), 



(26) 



which approximates the determinant of the continuous matrix 
operator defined over the Keldysh contour by a matrix deter- 
minant of the discrete matrix 1 + M. Hence, for each value 



of t that we wish to calculate the /-electron Green's function, 
we must first generate the corresponding matrix 1 + M for the 
Keldysh contour that runs out to time t and then take its deter- 
minant. Since the matrix can be generated solely from the pa- 
rameters U, t, T, and A(cj), this algorithm is easily parallelized. 
We first generate the function X(u>) on a discrete real frequency 
grid [by solving for the conduction electron Green's function, 
Eq. (|4}] on the master node, and then send that data to all of 
the slave nodes of the parallel process. Each slave process cal- 
culates the relevant matrix (which is a general complex matrix, 
with no special symmetries or properties), diagonalizes the ma- 
trix to find its eigenvalues, and then computes the determinant 
by taking a product of the eigenvalues. This is then sent back to 
the master who computes (t) from Eq. (II 01 and sends a new 
value of t to the slave to continue the computation. The algo- 
rithm has essentially a linear scale up in the parallelization, and 
it is quite efficient, since the limiting step is the diagonalization 
of the matrix, which is optimized by the local implementation 
of LAPACK and BLAS on the given parallel computer. We per- 
form the majority of our calculations on a Cray T3E, which lim- 
its the matrix sizes to approximately 2100 x 2100 on a 256 MEG 
node, and we can increase the matrix size slightly when work- 
ing on higher memory nodes, but the diagonalization time then 
becomes a limiter, if it takes longer than the queue limits for 
those nodes on a given machine. 

We fix the grid spacing on the imaginary axis, since we find 
the results are not too sensitive to changes of the step size there, 
and the value Atj ma g = 0.05 is sufficient for our purposes. 
On the real time axis, we generally take At rea j to vary from 
0.1 down to 0.0125. But since the calculations at different val- 
ues of t are completely independent of one another, we do not 
need to use the same grid spacing of the i-values for which we 
compute G^(t). We find that choosing a real time axis spac- 
ing of 0.2 or 0.1 [for generating the data for G^(i)] is usu- 
ally sufficiently accurate. We perform an Akuba shape pre- 
serving spline and evaluate the splined Green's function on a 
grid of size 0.01 or 0.005 before evaluating the cosine Fourier 
transform in Eq. ill 2b . This allows the accuracy to improve 
for larger frequency values, without producing a significant in- 
crease in computational time. Finally, the results of the Fourier 
transform, especially at small frequencies, depend on the cut- 
off or maximal time value where G^ (t) is evaluated. Since the 
Green's function decays at large times, imposing a cutoff is like 
replacing the Green's function by for times larger than the 
cutoff. We find that this usually causes no significant errors to 
the calculations when the maximum value of the Green's func- 
tion is less than approximately 10~ 4 to 10~ 5 in magnitude at 
the point where the cutoff is imposed. 

IV. Numerical Results 

When U = 0, the system is noninteracting, and the f- 
electron DOS is a delta function for all temperatures. When U 
increases, the /-electron DOS broadens into a regular function 
and picks up T-dependence (surprisingly, the conduction elec- 
tron DOS is always independent of temperature 1 1 1 ) . For small 
values of U, we expect the low-temperature /-electron Green's 
function to be a sharply peaked function with unit weight. The 
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Fig. 2 

Logarithm of the absolute value of (t) for U = 1 and T = 5. 
We plot results for At REAL = 0.1, 0.05, 0.025, and 0.0125. Note 
how the Green's function has a simple exponential decay at 

large times, and how the decay constant depends on the 
discretization size. the thick curves are the calculated data, 
and the thin curves are the extrapolated results. 



Fourier transform of this function to the real axis, will then be 
an exponentially decaying function, with a slower decay for a 
more sharply peaked function. Hence, we expect the greater 
Green's function to have an exponential tail for large times. At 
small times, the Green's function approaches p-t — 1 for all tem- 
peratures, and the curves for different temperatures start to sep- 
arate only for larger times. 

We illustrate the output of our calculations for the case U = 1 
and T = 5 in Fig.|2] We plot the logarithm of the absolute value 
of the real part of G/ (t) for four different choices of A£ rea j . We 
find, in all cases, that the Green's function has an exponentially 
decaying behavior at large times, so we append an exponential 
function out to large times, in order to have the Green's func- 
tion smaller than 10~ 4 at the maximal t cutoff. The extrapo- 
lated curves are represented by the thin lines. Note how there 
is a clear dependence of the Green's function on the step size 
taken along the real axis. Since this is a semi-log plot, it sug- 
gests that one extrapolate the logarithm of G^ (t) to determine 
the Ai rea j — > limit. But this procedure becomes problem- 
atic once the Green's function crosses zero, or has oscillatory 
behavior in the tails, so we do not carry out such a procedure 
here (note one might have expected (t) to have a quadratic 
dependence on Ai rea j due to a Trotter formula, but that does 
not apply here, since there is no simple Trotter breakup of the 
continuous matrix operator we are computing the determinant 
of). 

Next we examine the same set of parameters, but at lower 
temperature T = 0.15 in Fig. [3] Here the dependence on the 
discretization size is much stronger, with the exponential decay 
quite slow for the largest Ai rea j. This shows that the discretiza- 
tion size needs to be reduced as T is reduced, making lower 
temperature calculations much more difficult than higher tem- 
perature. Also, the maximal cutoff in time needs to be pushed 
farther out, unless one can append an extrapolated functional 
form (as we do here) for large times. 

The next step is to perform the Fourier transform (after splin- 
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Fig. 3 

Logarithm of the absolute value of (t) for U = 1 and 
T = 0.15. We plot results for At REAL = 0.1, 0.05, 0.025, and 
0.0125. Note how the Green's function has a simple 
exponential decay at large times, and how the decay constant 
depends strongly on the discretization size. the thick curves 
are the calculated data, and the thin curves are the 
extrapolated results . 



ing the time data) as in Eq. d!2i . We can then take the data for 
the DOS for different Ai rea j values and try to perform a point- 
wise (in u>) extrapolation down to Ai rea j = 0. Since we do 
not know how the curves depend on Ai rea j, we use an n-point 
Lagrange interpolation formula, which is a linear extrapolation 
for n = 2, a quadratic extrapolation for n = 3 and a cubic ex- 
trapolation for n — 4. By checking the different sum rules and 
the values of the Green's function at the Matsubara frequen- 
cies, we can examine the accuracy of different extrapolation 
schemes. Sometimes it is better to use all of the data and a 
large n Lagrange extrapolation scheme. Other times, the large 
step-size error is so big, that that data is not trustworthy, and 
it is more accurate to use an extrapolation with just the smaller 
discretization sizes (usually a linear extrapolation method with 
the smallest two values of Ai rea j is used). We call this extrap- 
olation technique the (^-extrapolation. 

Often, we find that the exact value for the Green's function 
at the lowest Matsubara frequency lies in between the value 
for the smallest Ai rea j and the value generated from the 6- 
extrapolation. In this case, we usually can improve the DOS if 
we perform a second extrapolation, averaging those two DOS 
to produce the correct value for Gf(iu>o). We call this extrapo- 
lation scheme Matsubara-extrapolation. 

To illustrate how the extrapolation procedures work, we first 
examine the high-temperature case T = 5. The results for the 
DOS for the four different discretization sizes and for the S- 
extrapolation is shown in Fig.|4] It is apparent from the figure 
that the DOS broadens and the peak is reduced as the discretiza- 
tion size is made smaller. The 6 extrapolation uses the two 
smallest values of Ai rea j and a linear extrapolation (we found 
that gave the most accurate results). A summary of the moment 
sum rules and the spectral formula for the Matsubara Green's 
functions is given in Table |ll] It is clear from that data that a 
systematic extrapolation is possible, and that the final result is 
highly accurate for the DOS (errors are less than 0.1% for the 
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DOS FOR U = 1 AND T = 5 AND DIFFERENT VALUES OF At REAL . ALSO 

PLOTTED IS THE RESULT OF THE ^-EXTRAPOLATION. ONE CAN SEE 
CLEARLY THAT THE DOS IS CONVERGING TO A WELL-DEFINED LIMIT AS 
THE DISCRETIZATION ERROR IS REDUCED. 



first moment, 0.03% for the second moment, and 0.003% for 
the lowest Matsubara frequency). It is hard to judge the ab- 
solute accuracy of the DOS from these integrated sum rules, 
but as a general rule, if the moment sum rules are accurate to 
better than 1% and the Matsubara frequency Green's functions 
are accurate to better than 0.1%, then the DOS has an absolute 
accuracy that is probably better than a few percent for most fre- 
quencies, except those near the tails, where the DOS is small 
and may even go negative. 

TABLE II 

Summary of moment sum rules and the Matsubara Green's 
function for the case u = 1 and t = 5. the first four rows are 
different values of at rea l- the last row is the exact result. 
Shift is the value added to the DOS to satisfy the sum rule in 

eq. int . 



Case 


Eq. G3 


Eq. CD 


G f {iw Q ) 


shift 


0.1 


-0.00787 


0.15805 


-0.063667 


0.00484 


0.05 


-0.00991 


0.19933 


-0.063660 


0.00500 


0.025 


-0.01087 


0.21863 


-0.063628 


0.00253 


0.0125 


-0.01137 


0.22873 


-0.063623 


0.00125 


(5-extrap. 


-0.01241 


0.25007 


-0.063600 


-10~ 6 


Exact 


-0.01240 


0.25 


-0.063598 


0.0 



We next examine the low-temperature case of U — 1 and 
T = 0.15 in Fig. |5] Note how the large At rea j case has 
an extremely narrow DOS, and how the DOS broadens dra- 
matically as the step size is reduced. We are able to make 
this reduction, since we can append the exponential tail out to 
large times, which allows us to perform the relevant Fourier 
transforms. We plot two extrapolation techniques: the 5- 
extrapolation, using a linear extrapolation with the two smallest 
step sizes (since that is the most accurate), and the Matsubara- 
extrapolation procedure, which averages the Ai rea j = 0.0125 
DOS with the 5-extrapolated DOS in order to properly repro- 
duce the lowest Matsubara frequency Green's function. Since 





12 


I 1 I 1 I 1 I 1 I 

1 At=0.1 




10 


f| At=0.05 

j'lu At=0.025 




8 


■| ■ At=0.0125 - 


*-h-> 


6 


6 — extrap 


y— 


4 


- ^^^^^^ 




2 












-0.4 -0.2 0.2 0.4 
Frequency a> [t*] 

Fig. 5 

dos for u = 1 and t = 0.15 and different values of at real . 
Also plotted is the result of the (5-extrapolation and the 
Matsubara-extrapolation. One can see that the DOS is 

CONVERGING TO A WELL-DEFINED LIMIT AS THE DISCRETIZATION ERROR 
IS REDUCED, BUT THE ACCURACY IS SIGNIFICANTLY REDUCED RELATIVE 
TO THE HIGHER-TEMPERATURE CASE IN FlG.|4] 



the (5-extrapolated result overshoots the correct answer, this av- 
eraging procedure significantly enhances the results. Even so, 
one can see from TablelHIIthat the errors are much larger at low 
temperature than at high temperature (this was already apparent 
from Fig.[3J. The error in the moment from Eq. dl5l is 1%, the 
moment from Eq. d!6i is 1 .5% and for the Matsubara frequency 
Green's function Gf(iuji) is 0.07%. 

TABLE III 

Summary of moment sum rules and the Matsubara Green's 
function for the case u = 1 and t = 0. 15. the first four rows 
are different values of at real . the last row is the exact 

RESULT. 



Case 


Eq. O 


Eq. d 


G f (iu> ) 


Gf(iu>i) 


0.1 


-0.05361 


0.12596 


-1.94404 


-0.6832 


0.05 


-0.08145 


0.20057 


-1.855512 


-0.6742 


0.025 


-0.09030 


0.21792 


-1.81501 


-0.6694 


0.0125 


-0.09760 


0.24410 


-1.79435 


-0.6668 


(5-extrap. 


-0.10307 


0.25936 


-1.77368 


-0.6642 


M-extrap. 


-0.10112 


0.25390 


-1.78105 


-0.6651 


Exact 


-0.10217 


0.25 


-1.78106 


-0.6656 



Hence these calculations become significantly more chal- 
lenging as the temperature is lowered. To understand the ther- 
mal evolution of the /-electron DOS, we summarize data col- 
lected for a number of different temperatures in Fig. [6] In addi- 
tion, we plot the temperature-independent conduction electron 
DOS with the magenta dashed line. Note how the conduction 
electron DOS is rather broad, but the /-electron DOS becomes 
sharply peaked at low temperature. By scaling the results to 
T = 0, we conjecture that the maximum of the /-electron DOS 
approaches 21 as T — > 0. As the DOS becomes more sharply 
peaked, we need to have a transfer of some spectral weight to 
larger energy as well (| ui | > 1.5), since the second moment sum 
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Fig. 6 

Summary plot for the /-electron DOS for U = 1 and a number 
of different temperatures. The magenta dashed line is the 

conduction electron dos (which is independent of 
temperature). Inset is a plot of 1/Af(u> = 0) versus T, to 

EXTRAPOLATE THE DOS DOWN TO T = 0. WE PREDICT THAT THE DOS 
GROWS TO A PEAK HEIGHT OF ABOUT 2 1 AT T = (BLUE TRIANGLE IN 

INSET). 



rule [in Eq. 1161 1 is independent of temperature, and as the peak 
grows in height, it contributes less to that sum rule. There is 
an interesting contrast in the DOS of the two different particles. 
The conduction electron DOS is broad and does not evolve with 
temperature, while the /-electron DOS has strong temperature 
dependence becoming sharply peaked at low temperature and 
weak coupling. By carefully performing extrapolations of the 
exponential decay of the greater Green's function at large time 
and of the discretization error of the DOS, we can produce ac- 
curate results for the /-electron DOS. The calculational needs 
can easily go beyond computational resources at low tempera- 
ture, though. 

Next we focus on the strong-coupling regime with {7 = 5. 
In this case, the conduction electron DOS has a correlation in- 
duced gap from the Mott transition, which occurs at U = 2. 
When the DOS develops a gap at low u>, the greater Green's 
function has a strong oscillatory component at large times. Un- 
fortunately there does not appear to be any simple rule that 
could be used to append an extrapolated result to the Green's 
function at large time. This greatly reduces the ability to calcu- 
late accurate results at low temperature, because the discretiza- 
tion size needs to be small, but the cutoff in time needs to be 
large, and this often goes beyond available computer resources. 

In Fig. [7J we plot the greater Green's function versus time 
for T = 1. In panel (a), the short-time results are shown — 
note how the curves for different discretization size lie almost 
on top of each other. In panel (b), we show the larger time re- 
sults (the smallest discretization size [green curve] only goes to 
t = 27). What is interesting to note is that the amplitude of 
the oscillations is reduced as the discretization size is reduced, 




20 25 30 
Time [1/t*] 

Fig. 7 

Greater Green's function for U = 5 and T = 1 for different 
At RE AL- Panel (a) is the short-time results and panel (b) is the 

LONGER-TIME RESULTS. NOTE HOW THE CURVES ARE VIRTUALLY 
IDENTICAL FOR SHORT TIMES, BUT THERE IS A SYSTEMATIC REDUCTION 
IN THE AMPLITUDE OF THE OSCILLATIONS AS At REAL IS REDUCED. 



but the period remains essentially the same. This might suggest 
to try to determine the amplitude reduction factor for the limit 
Ai rea j — > 0, and apply that to the curves for finite A£ rea j to ex- 
trapolate to the exact result. The problem is that such a scheme 
does not work well, as it tends to generate an unphysical peak 
in the low-frequency DOS, and it produces worse agreement for 
both the moments and the Matsubara Green's functions, so we 
don't discuss that scheme further. 

We are thus left with using the same 5 and Matsubara ex- 
trapolation schemes that we used for the small-[/ case. At high 
temperature, everything works out reasonably well, as can be 
seen in Table llVl But we do not achieve anywhere near as high 
an accuracy as for small U. When we go to lower temperatures, 
the accuracy becomes worse, and there is limited improvement 
from the extrapolation schemes. This occurs because when the 
discretization size is too large, the DOS becomes negative in the 
gap region, while when the size is small enough to correct the 
DOS in the gap region, it suffers from Gibb's oscillations due to 
the sharp cutoff in time, since the cutoff is not large enough for 
an accurate Fourier transform. These results are summarized in 
Fig.|8]and TablelVl Note that the sum rules do not change much 
as the DOS varies in the gap region, because the overall DOS is 
small and because we multiply by powers of u> which suppress 
the weight in the integral from the gap region. 

Finally, a summary of the DOS data for U = 5 is plotted 
in Fig. |9l We see that there is significant subgap DOS at high 
temperature, which is reduced as T is lowered. We also see the 
peak in the /-electron DOS move towards the band edge as T 
is lowered. Finally, it appears that the /-electron DOS and the 
conduction electron DOS will both share the same bandwidth 
at T = 0. We are severely limited by how low we can go in 



TABLE IV 

Summary of moment sum rules and the Matsubara Green's 
function for the case u = 5 and t = 5. the first two rows are 
different values of at real . the last row is the exact result. 
Shift is the value added to the DOS to satisfy the sum rule in 
EQ. (fl4l . 



Case 


Eq. O 


Eq. OSJl 


Gf(iw ) 


shift 


0.1 


-0.30191 


6.19393 


-0.062162 


0.00486 


0.05 


-0.30299 


6.21921 


-0.062147 


0.00502 


(5-extrap. 


-0.30396 


6.24237 


-0.062120 


0.00010 


Exact 


-0.30422 


6.25 


-0.062101 


0.0 
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Fig. 8 

DOS for U = 5 and T = 1. Three different values of At REAL 

AND THE 5 EXTRAPOLATION ARE SHOWN. IN THE LOWER PANEL, THE 
VERTICAL AXIS OF THE FIGURE IS BLOWN UP TO HIGHLIGHT THAT THE 
DOS IN THE GAP REGION IS NEGATIVE FOR LARGE DISCRETIZATIONS 
AND BECOMES POSITIVE ONLY AS THE DISCRETIZATION IS REDUCED. 



temperature and still maintain positive DOS in the gap region. 



V. Conclusions and Outlook 

In this contribution, we have presented a summary of numer- 
ical calculations employing a nonequilibrium formalism over 
the so-called Keldysh contour. The calculations involve deter- 
mining the determinant of a continuous matrix operator defined 
along the Keldysh contour, which is found by calculating the 
determinant of a discretized version of the operator, which is 
a general complex matrix, defined differently for each value of 
time. This formalism is easily parallelized, is efficient on each 
node, and has nearly linear scale-up. We examined the simplest 
problem with this numerical algorithm — the /-electron spec- 
trum of the Falicov-Kimball model. This problem is useful be- 
cause a number of sum rules exist, that allow one to determine 
the accuracy of the calculations. 

Our results show that one needs to carefully perform an anal- 
ysis of the dependence of the solutions on the discretization 



TABLE V 

Summary of moment sum rules and the Matsubara Green's 
function for the case u = 5 and t = 1. the first three rows 
are different values of at rea l- the last row is the exact 
result. Shift is the value added to the DOS to satisfy the sum 

RULE INEQ. fUl. 



Case 


Eq. {[3 


Eq. CD 


Gf(iu> ) 


shift 


0.1 


-1.01562 


6.28460 


-0.203556 


0.00504 


0.05 


-1.01334 


6.30267 


-0.203146 


0.00507 


0.025 


-1.01012 


6.29576 


-0.202731 


0.00504 


5-extrap. 


-1.00954 


6.30842 


-0.202735 


-0.00315 


Exact 


-1.00177 


6.25 


-0.203916 


0.0 
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Fig. 9 

Summary DOS for U = 5. The data for T = 5 comes from the 

^-EXTRAPOLATION, FOR T = 2 FROM At REA L = 0.05, AND THE LOWER 

TEMPERATURES HAVE At REA L = 0.025. WE ALSO INCLUDE THE 
CONDUCTION ELECTRON DOS WITH THE DASHED LINE. NOTE HOW THE 
SUBGAP DOS FOR THE /-ELECTRON DECREASES AS T IS LOWERED, AND 
HOW THE DOS MOVES MORE TOWARDS THE BAND EDGE AS T IS 
REDUCED. WE ALSO SEE THE LARGE |ti>| DOS START TO PINCH OFF AT 
THE CONDUCTION ELECTRON DOS BAND EDGE AS T IS LOWERED. WE 
EXPECT THE BANDWIDTHS OF BOTH THE CONDUCTION AND / ELECTRONS 
TO BE EQUAL AT T = 0. THE OSCILLATIONS APPARENT IN THE T = 0.8 
DATA ARE AN ARTIFACT OF A TIME-DOMAIN CUTOFF THAT IS TOO SHORT. 



to be generalized for each of these cases. One needs to find 
a way to self-consistently map an impurity or cluster problem 
onto the lattice problem (in the presence of the time-dependent 
field), and then perform similar computations along the Keldysh 
contour. The impurity-like problems are similar, but we need to 
perform a summation over the lattice wavevector (which can be 
represented by a double integral over a generalized joint den- 
sity of states for each matrix element of the contour-ordered 
Green's function) to complete the self-consistent algorithm. Far 
fewer sum rules will exist to benchmark the results, so an anal- 
ysis in terms of scaling with respect to the discretization size 
will need to be performed. One can check the nonequilibrium 
results in the linear-response regime with the results of Kubo- 
formula-based approaches, which will provide a stringent test 
of the quality of the numerics. 



along the Keldysh contour. Sometimes results can be extrapo- 
lated to the continuum limit, other times, this is not possible. 
We also find that these results often require a finer discretiza- 
tion at lower temperature. We anticipate similar numerical is- 
sues will arise in other nonequilibrium problems that employ 
the same formalism. 

These results provide useful benchmarks for more interesting 
nonequilibrium problems such as the interaction of a strongly 
correlated material with a strong external electromagnetic field 
or the nonlinear response of an inhomogeneous multilayered 
device to an external voltage (including a noise analysis of the 
current). These latter problems are likely to have more of an ap- 
plied interest to the military. The calculational formalism needs 
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